On the Effects of Different trans and cis Populations in Azobenzene Liquid Crystal Elastomers: A Monte Carlo Investigation

We investigate main-chain liquid crystal elastomers (LCEs) formed by photoresponsive azobenzene units with different populations of trans and cis conformers (from fully trans to fully cis). We study their macroscopic properties as well as their molecular organization using extensive Monte Carlo simulations of a simple coarse-grained model where the trans and cis conformers are represented by soft-core biaxial Gay–Berne particles with size and interaction energy parameters obtained by fitting a bare bone azobenzene moiety represented at atomistic level. We find that increasing the fraction of cis conformers, as could be obtained by near-UV irradiation, shifts the nematic–isotropic transition to a lower temperature, consistently with experiment, while generating internal stress in a clamped sample. An analysis of pair distributions shows that the immediate surroundings of a bent cis molecule are slightly less dense and more orientationally disordered in comparison with that of a trans conformer. Comparing nematic and smectic LCEs, actuation in the smectic phase proved less effective, disrupting the smectic layers to some extent but preserving orientational order of the azobenzene moieties.


■ INTRODUCTION
Stimuli-responsive materials are an important class of smart materials capable of responding to an externally applied trigger through an observable property change. These triggers can include temperature variation, application of an external electric or magnetic field, and irradiation with light, of interest here. 1 The latter applies, for example, to photoresponsive polymers 2,3 containing azobenzene moieties whose molecular shape anisotropy changes significantly upon irradiation with near-ultraviolet (UV) light causing a trans to cis isomerization. 4,5 If these moieties are themselves mesogenic and exhibit nematic or smectic phases, 6,7 it is found that their orientational, liquid-crystalline ordering can be decreased or destroyed upon irradiation with light even at constant temperature (with an effect similar to that of increasing temperature). 8−14 Somewhat amazingly the conformational changes take place even in very viscous glassy materials, possibly through photofluidization, 15 although the issue is still debated. 16 Photoresponsive liquid crystal elastomers (LCEs), characterized by a pronounced strain-alignment coupling, can be obtained when azobenzenes are covalently linked to a cross-linked polymeric network 1,17−19 or even just dissolved as guest in the LCE matrix. 20 Such materials can exhibit significant elastic deformations upon phase transitions between the orientationally ordered liquid crystalline state and the disordered isotropic phase and can thus act as photoactuators converting light energy into mechanical work. 20−23 In view of this, photoresponsive LCEs are promising candidates for a variety of applications, 24 including sensors and actuators, 1,25−27 photomechanical devices, 19 morphing surfaces, 22 soft robotics, 28,29 and more, also thanks to advances in chemical synthesis. 30 The photophysics process of switching between the trans and cis states (see Figure 1) of an individual azobenzene molecule has been modeled in many studies in vacuum 4,5 or, at the atomistic level, for azobenzene in isotropic solvents, 31,32 but simulations performed at this level of detail can be considered only for systems with a small number of molecules. Atomistic and multiscale modeling has also been performed for the low molar mass azobenzene mesogen 4,4′-dioctyloxyazobenzene (8AB8) to investigate its phase changes. 33 7AB (p,p′diheptylazobenzene) in 8CB (4-octyl-4′-cyanobiphenyl) has also been studied, in particular to examine the preferred location of trans-and cis-7AB in the smectic A phase of 8CB and finding that while the trans form is hosted inside the smectic layers, the cis tends to reside between the layers. 34 Recently, atomistic simulations of azobenzenes in molecularand polymer-network glasses have been carried out to study their behavior during and after photoactivation in a rigid matrix environment. 35,36 As far as liquid crystal polymers are concerned, apart from theoretical efforts 37 aimed at understanding azobenzene-based polymers, a few particle-based computer simulation studies have also been performed at the coarse-grained level, 38,39 particularly with the purpose of studying the photoselective azobenzene excitation followed by effective reorientation of the director under linearly polarized UV light excitation. 14 It is worth noting that, although very interesting in themselves, the details of the trans−cis conversion are not of central interest in terms of properties of the material exposed to isotropically distributed unpolarized UV light, if only for a time-scale argument. Indeed, the photoinduced trans-to-cis conversion is a very fast one that can even be performed in a few femtoseconds with a sufficiently powerful light source, while the rate of the spontaneous thermal cis-to-trans back conversion depends greatly on the specific compound 40 and, unless pumped back by irradiation in the visible, where the typical cis absorption peak occurs, is typically slow to very slow (even well over a day 41 ). We can then imagine that a system, initially with all azobenzenes in the trans state, could be suitably irradiated to yield a system with a certain controlled fraction of cis and that this system will be sufficiently long-lived to establish different molecular organizations according to the achieved fraction. Accordingly, our aim here is to investigate how the presence of this mixture of conformers will affect LCE phase transitions, sample deformation, and the differences of local structure around the two types of conformers. As far as we know, these aspects have not been dealt with before.
In the next sections, we describe our coarse-grained model for azobenzene LCEs and the range of Monte Carlo (MC) simulations performed on LCEs with various percentages of trans and cis conformers, monitoring ordering and structure as a function of temperature. Stress−strain curves at selected temperatures and trans and cis conformer ratio will also be provided. Conclusions will be collected in the last section.

■ MODEL AND SIMULATIONS
We choose a generic particle model where only the essentials of an azobenzene-based LCE are retained. This type of approach is probably unavoidable, considering that the actuation in LCE is a collective phenomenon involving a large number of particles. Moreover, reliable microscopic modeling of any cross-linked LCE system should involve samples whose size well exceeds the characteristic distance between the cross-links (that, in turn, significantly exceeds the size of a single molecule) to exclude the effects related to local irregularities in the LCE polymer network. Currently, sufficiently large sample sizes can only be achieved at the coarse-grained simulation level. Ideally, one would want to connect the benefits of the atomistic and coarse-grained treatment in a suitable multiscale approach. For this reason, in this study, we start from an atomistic-level optimization of the molecular structure for both azobenzene conformers, trans and cis, and use them to estimate the corresponding interaction parameters of the soft-core Gay−Berne (GB) attractive− repulsive potential acting between biaxial ellipsoidal particles representing the molecules 42,43 that is commonly used in largescale simulations of liquid-crystalline systems. 7 With the thus obtained GB molecular building blocks, we assemble suitable LCE networks, following the methodologies developed in refs 44−46 and perform large-scale MC simulations. From the results, we determine the shift of the nematic−isotropic (NI) transition temperature upon increasing the cis particles content, together with the magnitude of the internal stress generated upon irradiation of the sample with UV light. Moreover, we examine the positional and orientational intermolecular correlations in the vicinity of both trans and cis particles.
We have previously modeled and simulated LCEs using soft GB units 47 linked together by finitely extendable nonlinear elastic (FENE) 48,49 spring-like bonds. 44−46 The soft-core modification of the original GB (Appendix A) 42,50 reduces its harsh repulsion, making entangling of the chains more manageable, and turns out to be essential for a proper The dashed circle indicates one of the trans or cis azobenzene units whose united-atom representation is shown together with their coarse-grained GB ellipsoids in (b), trans, and (c), cis. The grey and blue spheres represent the −CH− groups and the nitrogen atoms, respectively. The conformational change at the double bond between the nitrogen atoms results in a significant modification of the fitted coarse-grained molecular shape depicted by the transparent ellipsoids. Insets: locations of bond site positions at the GB ellipsoid surface and easy bond directions for the main-chain bonding (red) and for the optional cross-linking via the equatorial site (blue), as used in simulation. For the trans isomer, the head and tail bonding sites are assumed to impose an easy direction along the molecular long axis, while the optional site on the equator is taken to impose an easy direction along the second-longest molecular axis. In the cis isomer, the linking directions of the head/tail sites are tilted to make an angle of ≈46°with the long molecular axis. equilibration of cross-linked GB polymer networks. In addition, a soft biaxial version of the GB potential (SBGB) has also been developed. 43 It is thus natural to extend this approach to photosensitive LCEs by replacing some or all of monomer units with SBGBs modeling either trans or cis azobenzene moieties (see Figure 1a).
It is worth noting that an atomistic model, apart from being hardly feasible in terms of resources, as already mentioned, could even be not ideal, since it has to pay the appealing feature of being chemically realistic with the specificity of being applicable only to a certain chemical formulation. Since the actuation phenomenon we aim to better understand takes place for a variety of molecular structures, which have in common only the presence of the photoresponsive azobenzene moiety, it is desirable to develop a minimal model where only the trans−cis conversion is accounted for and other chemical details abandoned.
In practice, we study, using MC simulations, the effect of replacing a certain fraction of trans azobenzenes with cis conformers in the LCE. We shall go from full trans to full cis, considering fractions of cis w c = 0, 0.2, 0.4, 0.6, 0.8, and 1.0 fixed during the simulation at various temperatures. In order to maximize the order-disruption effect and, consequently, the actuation, here all GB ellipsoids, both those of the polymer strands and the cross-linkers, are taken to be switchable azobenzenes as, for example, in ref 21. The coarse-graining step from the atomistic to coarse-grained description of a model azobenzene molecule is sketched in Figure 1b,c. As we can see, our minimal molecular model of an azobenzene photoresponsive unit consists of two aromatic rings connected via a N�N double bond whose loosening under the effect of UV light gives rise to the trans−cis isomerism 5,51 that is responsible for the difference in the effective molecular shape. At the coarse-grained level, both isomers are represented by biaxial ellipsoids (since no molecule is perfectly uniaxial), whose interaction parameters need to be determined from an atomistic representation. Since we now have various types of units, interactions between dissimilar biaxial GB particles have to be considered. 52 The main input parameters of the soft-core GB potential describing the two isomers are particle dimensions expressed through ellipsoidal axes (σ x , σ y , and σ z ) and the energy parameters (ϵ x , ϵ y , and ϵ z ) related to the attractive potential well depths for corresponding homogeneous particle alignments. 50 Their values are reported in Table 1 and have been obtained from atomistic representations at united atom (UA), rather than at all-atom level. Note that not only the molecular aspect ratio and the relative characteristic potential well depths can be derived from the fitting but also their absolute values that eventually translate into macroscopically observable quantities. More details can be found in Appendix A and Appendix B.
We consider a main-chain system where the elongated mesogenic moieties belong to the polymer backbone. In order to build an elastomer network, the fitted SBGB particles are bonded with other ellipsoids to build the main chains, as well as connected to other chains via cross-linking particles. In our minimal model, we neglect, apart from the change of shape, any realistic chemical details and just assume that each trans particle is decorated with a "head" and "tail" bonding site, that will be displaced at an angle in the cis molecule, following the conformational change. We also arbitrarily assume (ignoring at this idealistic model level any chemical difficulty) that the cross-linking particles have an additional bonding site on the equator. The positions of the bonding sites and the easy (preferred) directions of the resulting bonds, as assumed in the GB simulation, are indicated in the insets of Figure 1b,c.
The bonds themselves are modeled using the FENE potential 48 applied both to bond stretching and bending, 44 see Appendix A. The FENE potential, which is often used in modeling polymer chains, has the important feature of being stiffer for large deformations than its Hookean counterpart, imposing a maximum stretching and bending range.
All samples considered in this study contain N = 64,000 GB particles. Like in our previous work, 44−46 we consider LCE swollen with a liquid-crystalline solvent (in practice a monomeric excess), which are characterized by a greater mobility of polymeric chains and have been experimentally studied. 53−55 Thus, only one-half of the mesogens is used to build the elastomeric polymer network itself, while the rest represents the nonbonded swelling monomers introduced to facilitate the relaxation of polymer chains in simulation, while conveying nematicity at the same time. 44 Periodic boundary conditions 56,57 are applied in each sample preparation and in simulations to mimic a larger bulk sample. The polymer network is assembled following the steps similar to those presented in ref 46 to build polydomain LCE samples. Here, however, again to maximize actuation, we consider monodomain samples where the average molecular orientation, the nematic director, is the same throughout the sample and is imprinted parallel to, say, the laboratory z-axis. This is accomplished by growing, at low density, straight individual chains of particles, polydisperse and on average ∼12 monomers long, each starting from a random position within a cubic simulation box, directed at a random angle smaller than arccos(0.6) ≈ 53°, measured from the z-axis. Then, the active heads and tails of each chain are connected laterally to a nearest equatorial site on another, neighboring chain. A schematic representation of the typical resulting elastomer architecture is shown in Figure 1a. Finally, the system is soaked with swelling monomers and isotropically compressed almost to close-packing, increasing the density more than 20 times, to yield a cubic sample that we consider as a "reference" one. In principle, one should consider several of these reference samples to obtain a collection of LCE networks grown under the same conditions, but this is practically unfeasible in view of the computer time required. We have thus considered a single reference, also because, based on our previous experience, 44 −46 one replica for each cis particle concentration is sufficiently representative for a sample size as large as our one (N = 64,000) where the simulation box side significantly exceeds the average distance between the network cross-links.
Note that in our simulations, the reduced number densities where V is the sample volume and σ 0 is the shortest axis of the trans ellipsoid, depend on the cis particle fraction w c . This is because the fitted coarse-grained molecular volume (proportional to σ x σ y σ z ) is somewhat smaller for the cis isomer than for the trans one (see Table 1), and it is necessary to avoid shifts of the NI transition due to a mere dilution effect. Consequently, the final sample densities increase with increasing w c . In a simplified approach, we determine these densities, summarized in Table 2, by keeping the molecular packing fraction unchanged for all w c . Note that the ρ* values employed are taken to be relaxed equilibrium values that are assumed constant during a simulation run. At variance with this, experiments 22 have shown a reduction in sample density upon sample irradiation, which, however, is a transient phenomenon beyond the scope of the present study. To find equilibrium configurations at a given w c and temperature T, we performed MC simulations following the standard Metropolis algorithm. 7,58 In each MC cycle (a cycle being a set of N attempted MC moves), particles are selected one by one in random order and undergo an individual trial move that is accepted if the associated interaction energy change ΔU is negative, or with a probability of exp(−ΔU/k B T) otherwise, where k B denotes the Boltzmann constant. For a full definition of the interaction energy U, see Appendix A.
In each MC cycle, different trial move types are attempted separately: translational moves, rotational moves, 59 and bonded pair rotations. 44 In addition, every 5 MC cycles, a collective resize move is performed, where the cuboidal simulation box itself is randomly and affinely deformed at a constant volume, and subjected to the same Metropolis acceptance criterion as above (more details can be found in refs 44−46). The trial move amplitudes are adjusted on the fly to maintain a trial move acceptance ratio of around 50% while ensuring that the system is evolving in terms of particle positions and orientations. In the calculations of interaction energies, the computational effort is reduced by setting up particle neighbor lists 60 with interaction and list cutoffs at 6σ 0 and 7σ 0 , respectively. Moreover, neighbor list updates are facilitated by means of cell-linked lists, 56 and further parallelization schemes are employed to provide additional simulation speedup.
As a result, several 10 7 MC cycles have typically been performed for sample equilibration and ∼10 6 cycles for production. Typical averages of interest include (i) the reduced sample length λ z along the z-axis, the nematic director of the monodomain, to monitor actuation, (ii) the standard nematic order parameter P 2 to follow the underlying changes in orientational order, and (iii) the specific heat of the sample to facilitate the detection of the associated phase transitions. (i) Here, λ z is conveniently defined as the average simulation box length along the z-axis, divided by the side of a reference sample containing trans particles only. (ii) The nematic order parameter calculation is performed calculating, for a certain MC cycle, the nematic ordering matrix where u i stands for the long molecular axis of the ith molecule and I is a 3 × 3 identity matrix, and diagonalizing Q to identify its largest eigenvalue, which is then averaged over MC cycles to yield the order parameter indicated here as P 2 . 7 (iii) Further, the reduced specific heat (per particle and in units of k B ) is calculated as c* = d⟨U*⟩/dT*, where the reduced internal energy and temperature are introduced as U* = U/Nϵ 0 and T* = k B T/ϵ 0 , respectively, and ⟨···⟩ stands for averaging over MC cycles. (Here, ϵ 0 is the energy parameter related to the attractive GB potential well depth for a pair of perfectly aligned biaxial trans particles approaching along their shortest axes. 50,52 ) Finally, to identify layered smectic phases, the structure factor S(q) can be calculated as 61 where q denotes the scattering vector and r j stands for the center-of-mass position of the jth anisotropic bead. Then, to assess the degree of smectic ordering,

■ RESULTS AND DISCUSSION
In an aligned nematic containing trans azobenzene moieties either as part of the mesogens or as solutes, the presence of cis conformers represents a source of disruption in the alignment of their environment. This destabilizes the nematic phase and decreases the NI transition temperature, as observed experimentally. 8−14 Therefore, we first explore the phase diagram of the LCE samples we have generated, as a function of reduced temperature T* and cis particle content, w c , even for very large values of the latter. Simulation runs at all temperatures have been launched from the corresponding "reference" samples, except for a few cases at low temperatures (T* ≤ 3.0) where the runs were started from a pre-equilibrated configuration at a nematic temperature to facilitate MC equilibration. In an all-trans sample (w c = 0), plots of the sample dimension λ z and of the nematic order parameter P 2 against temperature ( Figure 2) reveal an isotropic-to-nematic phase transition at * T 7.5 NI . (Taking ϵ 0 ≈ 1.3 × 10 −21 J, as obtained from the UA fitting, this yields a NI transition at ∼430°C.) Below this temperature, the sample acquires an elongation along the z-axis, as imprinted upon sample preparation, that coincides with the nematic director, and a nonzero degree of nematic order (Figure 2b,c). The transition itself is also accompanied by a peak in the temperature dependence of the reduced specific heat, c* (Figure 2a). For this all-trans sample, there is also an additional peak in c* below T* ≈ 1.5 corresponding to the phase transition from the aligned nematic into the aligned and layered smectic phase�characterized by a nonzero smectic order parameter τ, Figure 2d�suffering from extremely slow MC relaxation.
As the cis particle content is increased, parallel nematic alignment is destabilized and the NI transition shifts toward Note, however, that a comparison with experiments is not straightforward here since both the NI transition temperature and its shift depend not only on the precise chemical nature of the mesogen but also on other features of the elastomer network and are thus quite system-dependent.) In a simulation, it is also possible to investigate the all-cis limit (w c = 1.0), which may be rather far from reach with photoirradiation experiments, and from a nonzero P 2 value (see Figure 2c), one can deduce that an aligned phase with orientational ordering is established even in this case. (Our simulations suggest that this happens for T* < 3.5 also in a nonbonded liquid-crystalline system consisting of monomers alone. In an LCE sample, such order is further enhanced by the mechanical aligning field 18 imprinted along the z-axis upon sample preparation. ) We note that while working with a certain fraction w c of cisazobenzenes, the specific selection of cis particles is chosen (randomly) before the simulation and then left fixed. Ideally, one could consider many sample replicas with different cis− trans particle sets and perform an average over simulations of these isomeric replicas. In practice, each simulation for our sample size takes too long to make this feasible. However, to check the robustness of our results, typically obtained from a single isomeric replica, we have repeated the temperature-scan experiment on a second alternative sample with a different cisparticle distribution (and the same w c ) and the results for the two replicas considered (especially the * T NI values) turn out to be in very good agreement (see Figure 2 and compare full and empty symbols).
The presence of shorter cis particles locally disrupts the nematic orientational order established by the more anisotropic trans particles. To gain more insight into how this happens, it is instructive to study the surroundings of each isomer type separately by calculating their radial distribution functions g(r ij ) and orientational order correlation functions G 2 (r ij ), for example, from a single molecular snapshot. g(r ij ) is calculated by identifying and counting all particle pairs where the first particle belongs to the selected azo-species and then creating a corresponding histogram with respect to the interparticle distance r ij . The histograms are normalized with respect to an uncorrelated uniform fluid and then averaged over all possible selections of the first particle of the given species. As such, g(r ij ) is then proportional to the local density of all particles surrounding a particle of a given azo-species. A similar analysis can be performed by averaging P 2 (u i ·u j ) over eligible pairs of particles i and j (instead of just counting the pairs), where u i and u j denote the particle long axes and P 2 (x) = (3x 2 − 1)/2 is the second Legendre polynomial, giving the second-rank orientational pair correlation 7 The correlation properties in a sample with 20% cis-content are summarized in Figure 3a. The radial distribution functions g c (r ij ) and g t (r ij ) for the surroundings of cis and trans isomers (top plate), respectively, show no major qualitative differences at a nematic temperature of T* = 3.0, but the details are compatible with the somewhat different particle shapes (see Table 1). The difference of average local particle densities in the vicinity of cis and trans particles, given by Δg(r ij )ρ*, where Δg(r ij ) = g c (r ij ) − g t (r ij ) and ρ* is the average sample density, is shown as inset in Figure 3a at different T* in the aligned nematic phase. At smallest distances, g c (r ij ) > g t (r ij ) due to the smaller σ x for the cis isomer, which is followed by a pronounced g c (r ij ) < g t (r ij ) region. As a result, within a small sphere ∼4σ 0 in diameter surrounding a cis particle, there are slightly less molecules than in an equal-sized sphere surrounding a trans particle. This appears to be consistent with observations and hypotheses by Liu. 63 Recall that here w c = 0.2, so that the vast majority of isomers are in trans state. A cis particle thus represents a slight ordering disruption and is somewhat less efficiently embedded into its predominantly trans-environment, than any trans particle itself. From the inset of Figure 3a, it also follows that these effects are more pronounced at lower temperatures.
Further, the orientational correlation functions G 2 (r ij ) are shown for w c = 0.2 in the bottom plate of Figure 3a separately for both isomers. At small r ij , they show a similar structure as their g(r ij ) counterparts, while for large r ij where the shortrange correlations are already lost, they are expected to relax to

ACS Applied Polymer Materials
pubs.acs.org/acsapm Article a value equal to the product of the nematic order parameters of the central (cis or trans) and a far-away (any type) azobenzene particle. (This follows from the addition theorem for secondorder spherical harmonics in the absence of biaxial molecular ordering.) From the lower G 2 (r ij → ∞) value for a cis particle at the center in comparison with the corresponding value for a central trans particle, one can conclude that the degree of orientational order of cis particles themselves (in a trans-rich w c = 0.2 sample) is lower than that of the trans particles. Similar behavior is also systematically observed at larger w c values (see Figure 3b): the longer trans particles are always more orientationally ordered than the shorter cis ones, which results in a higher G 2 (r ij → ∞) value in case of a central trans particle, compared with that of a central cis particle. At the same time, the G 2 (r ij → ∞) values, at given temperature, for both particle types decrease with increasing w c due to a general weakening of the nematic ordering (Figure 2c). Studying photoresponsive LCE, it is also interesting to look at sample behavior when increasing the w c fraction at fixed temperature. Figure 4a shows a series of molecular snapshots in such a situation at T* = 4.0 where the actuative effect upon changing w c is largest (as can be deduced from Figure 2). At this temperature, an all-trans (w c = 0) sample is nematically ordered and the sample strongly elongated along the z-axis. The increase of cis-content, w c , results in a gradual destruction of orientational order, as well as in the loss of sample elongation. This behavior could be exploited in actuation to extract useful mechanical work from the system. The simulation box side λ z vs w c is shown as inset, Figure 4b, for different temperatures. All dependences decrease monotonically, whereby lower temperatures generally result in a higher value of λ z .
At low enough temperatures, the simulated LCE samples are in the smectic phase. Given that in an LCE system, the ellipsoidal molecules are bonded in a rather complex crosslinked network, the smectic layers are well-defined but present plenty of defects. The molecular snapshot close-ups depicting actuation as w c increases in a simulated smectic LCE at T* = 1.0 are shown in Figure 5a. The increasing presence of the short cis particles additionally perturbs the smectic layers initially consisting of longer trans particles, but the orientational order, together with sample elongation, persists even in an all-cis sample at w c = 1.0. In addition to molecular organization snapshots (Figure 5a), the presence of positional smectic ordering is further evidenced by the structure factor patterns depicted in Figure 5b where a series of peaks along the z-axis (parallel to the smectic layer normal and the nematic director) is clearly visible. (The less pronounced arcs in the direction perpendicular to z are a signature of interlayer molecular ordering.) Figure 5c further shows the evolution of the smectic order parameter τ with increasing cis-content: at first, τ decreases due to the somewhat disrupted layering. However, approaching the all-cis sample, τ increases again, which can be attributed to the absence of the perturbing longer trans particles in a cis-rich environment. Note also that with increasing w c , the structure factor peaks move to a slightly higher wavevector value |q 1 |, compatibly with the shorter long axis of cis particles.
Above we have seen that upon light-induced actuation, an initially elongated sample contracts along the nematic director, z-axis. Such behavior applies to a free (unclamped) sample. Alternatively, one could also consider an irradiation experiment with a clamped sample whose length λ z is fixed. A trans− cis photoisomerization under such conditions is inevitably accompanied by a generation of internal stress resulting from an increase of molecular disorder that is, in turn, balanced by the external stress of the clamps. To estimate the magnitude of stress thus generated at a given temperature, we consider a contracted sample with w c ≠ 0 and apply external stress to stretch the sample along the z-axis to its initial length λ z measured at w c = 0. Each stretching experiment is performed through an isostress MC run 57 where, apart from temperature, engineering stress Σ zz along the z-axis (here calculated with respect to the surface area corresponding to a face of the alltrans reference sample) is fixed. Having fixed an additional extensive thermodynamic variable, the Metropolis acceptance

ACS Applied Polymer Materials pubs.acs.org/acsapm
Article criterion has to be modified: the internal energy change ΔU is replaced by the change of enthalpy ΔH, 57,64 where the interaction enthalpy H is given by U − Σ zz Vλ z ; V here denotes the corresponding sample volume. The reduced external engineering stress is conveniently defined as * = / zz 0 3 0 . It is quite impossible to guess the correct value of stress σ* necessary to stretch an irradiated sample far enough to exactly match its initial length. For this reason, we simulate full stress− strain curves λ z (σ*) and perform the necessary extrapolation between the data points thus obtained to find σ*. Each of the necessary simulation runs was started from a sample previously equilibrated at the required conditions (T*, w c ) and zero stress. In Figure 6, stress−strain curves for two values of w c (0.2 and 0.4) and two temperatures (4.0 and 6.0) are plotted. In the low-stress region, σ* < 0.03, the dependence is close to linear and, as a byproduct, the corresponding isothermal Young elastic moduli can be estimated from the λ z (σ*) slopes The modulus decreases with increasing cis particle content and decreases with temperature; its values are relatively high but compatible with our previous results for simulated LCE in the nematic phase.
Turning now to the values of internal stress σ* in a clamped sample, at T* = 4.0, we obtain 0.010 and 0.017 for w c = 0.2 and w c = 0.4, respectively. At T* = 6.0, the corresponding values are 0.009 and 0.021. (Assuming the same values for ϵ 0 and σ 0 as above, the given σ* values are to be multiplied by ∼30 MPa to obtain the engineering stress Σ zz .) At a given temperature, internal stress in a clamped sample thus increases with irradiation (leading to higher cis-content w c ) since a larger sample contraction has to be compensated for by the clamps.  Smectic layering is partially disrupted as the cis-particle content increases, but the orientational ordering and sample elongation are preserved even for w c = 1.0. In an all-cis sample, the smectic ordering is to some extent restored.

ACS Applied Polymer Materials pubs.acs.org/acsapm Article
For the same reason, at a given w c , one would expect internal stress to increase also with increasing temperature, but this is not always the case since at higher temperature our samples also become softer, which reflects in a reduction of the Young's elastic modulus. Finally, the typical magnitude of the generated internal stress is estimated to be of the order of several hundred kPa.

■ CONCLUSIONS
We have performed a computational study of the effects that different percentages of trans and cis azobenzene monomeric units can have on the properties of a model LCE. More specifically, a large number of coarse-grained MC computer simulations (over 80) have been employed to study shape changes in monodomain main-chain LCEs formed by azobenzene moieties with various fractions of trans and cis conformers. The photoresponsive moieties were represented by soft-core GB biaxial ellipsoids whose interaction parameters were obtained from a fit of the UA-level interaction energy for selected pair configurations. Our results reveal a decrease of the NI transition temperature upon increasing cis population (e.g., by increased sample irradiation). In samples with prevailing trans-particle presence, the immediate neighborhood of cis particles is slightly less dense than that of the trans ones that are more compatible with their average environment.
Simulations of light-induced actuation at fixed temperature reveal a significant contraction of the sample along the monodomain director, which could be exploited to perform useful mechanical work. This effect is most pronounced when actuating at temperatures not too far below the NI transition. Finally, from the simulated stress−strain curves, we conclude that in a clamped sample, internal mechanical stress of the order of several hundred kPa can be generated upon irradiation with light and that this stress increases with the content of cisisomers, but not necessarily with temperature.
■ APPENDIX A

Coarse-Grained Pair Potentials
As mentioned in the text, the coarse-grained pair potentials considered in our study are the soft-core biaxial GB (SBGB) for the nonbonded interactions and the finitely extensible nonlinear potential (FENE) for the bonded ones.
The biaxial soft-core GB potential 43 acting between dissimilar ellipsoidal particles i and j is defined as 47 is the linear repulsive soft-core potential and U ij GB is the biaxial GB potential. 50,52 The switches smoothly between the soft-core repulsion at small interparticle distances r ij and the attractive potential well of the GB potential at larger r ij . Above, m and n are parameters (we use m = −70ϵ 0 /σ 0 and n = −100/σ 0 44−46 ), and σ ij is the anisotropic distance function for the two particles given by = r A r (2 ) ij ij where r ij is the interparticle vector and r ij = |r ij |. The distance function σ ij also depends on the relative particle orientation, which is encoded into the A i j overlap matrix defined by , where M i is the rotation matrix that brings any vector from the laboratory reference frame to that of molecule i. 65 Moreover, S i are range matrices describing the molecular shape, with the corresponding σ x , σ y , and σ z ellipsoidal axes on the diagonal and zeros elsewhere.
The biaxial GB potential is now given by where σ ij and r ij were defined above and σ c denotes the minimum interparticle contact distance. (In the case of dissimilar particles, an arithmetic average of the shortest axes of the particles is taken. 52  We employ the FENE potential 48 for interparticle bonds to describe both their stretching and bending. Consider a bond connecting two neighboring ellipsoids i and j, with a length of s ij . Then, the deviation from its equilibrium value s e equals δs ij = |s ij − s e |, and let us assume that δs ij < δs m . (There is a limitation on stretching.) Similarly, the bond angle θ ij is defined as the angle between the preferred bond directions imposed by each of the particles. 66 The deviation from its equilibrium value θ e is δθ ij = |θ ij − θ e |, with a requirement that δθ ij < δθ m . (There is a limitation on bending, too.) Now, the FENE bond energy is given by 66  We performed a simple energy minimization of atomistic structures of the two isomers of azobenzene with "Argus Lab" 68 choosing the semiempirical Quantum Chemistry (QC) method AM1 (Austin model 1) and the Restricted Hartree− Fock self-consistent field method (for closed shells). 69 We recall that we are only interested in getting the basic shape change upon trans−cis photoisomerization of the bare bone azobenzene unit in Figure 1, thus using more sophisticated QC methods is unnecessary, if not inappropriate. Using the force field published in ref 31, we have obtained molecular ellipsoids overdimensioned with respect to the atomistic structures. Hence, we have applied a UA approach, maintaining the same nomenclature of the force field but including into the carbon atom of the benzene ring the hydrogen attached to it. 7 Atomic masses and charges have been modified accordingly to take into account the collapsed hydrogens. The dimensions of the biaxial molecular ellipsoid are computed using the "shape tensor" T s , a modified version of the inertia matrix in which the masses are substituted with the van der Waals radii 70 to consider steric hindrance of the atomic centers, as reported below i k j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z Here, r i is the van der Waals radius of the ith atom and the total sum of radii is denoted by r T = ∑ i r i . The eigenvalues λ xx , λ yy , and λ zz of T s are related to the dimensions of the ellipsoid and can be used to define its axes as follows When the molecule is flat (as in the trans isomer), the length of the shortest axis approaches zero. Therefore, the standard distance of π−π stacking, 0.35 nm, has been used instead.
The energy parameters ϵ x , ϵ y , and ϵ z related to the potential well depths are computed using the atomistic structure and the UA force field described above. Two identical molecules have been placed at different distances, and their interaction energy within the chosen force field has been computed. The interaction energy has been obtained as a function of intermolecular distance for the three most significant pair configurations: side-by-side, face-to-face, and end-to-end. These profiles have then been fitted with the soft-core GB potential, using the ellipsoid dimensions obtained from the shape tensor.
The parameters summing up fitting results for the cis and trans UA azobenzene conformers represented as prolate GB ellipsoids are reported in the text (